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A  Three-Dimensional  Finite-Difference  Time-Domain  Formulation 
For  Nonlinear  Materials  with  the  Frequency-Dependent  Electric  Conductivity  and  Polarization 

S.  Joe  Yakura  and  David  Dietz 
Air  Force  Research  Laboratory,  Directed  Energy  Directorate 
Kirtland  AFB,  NM  871 17-5776 

Abstract 

In  this  paper,  we  present  a  three-dimensional  finite-difference  time-domain  ( FDTD)  algorithm  that  is  used  to 
evaluate  the  propagation  of  electromagnetic  waves  in  conductive  and  dispersive  materials  that  exhibit  the 
frequency-dependent  electric  conductivity  and  polarization.  We  consider  a  case  where  the  electric  conductivity  has 
the  linear  property ,  specifically  through  the  first-order  (linear)  electric  conductivity  function ,  and  the  electric 
polarization  has  both  linear  and  nonlinear  properties ,  specifically  through  the  first-order  (linear)  and  third-order 
(nonlinear)  electric  susceptibility  functions .  The  resulting  FDTD  algorithm  shows  that  the  nonlinear  dispersive 
material  with  a  third-order  susceptibility  function  results  in  coupled  nonlinear  cubic  equations  for  the  three 
components  of  the  electric  field  vector ;  relating  the  next-time-step  electric  field  vector  to  the  previous-time- step 
electric  field  vector.  This  contrasts  the  usual  algorithm  of  the  linear  conductive  and  dispersive  material,  which  has  a 
simple  linear  relationship  between  the  next-time- step  electric  field  and  the  previous-time-step  electric  field. 
Consequently,  the  coupled  nonlinear  cubic  equations  must  be  solved  at  each  time  step  to  simulate  the  behavior  of 
the  electric  field  vector  in  the  nonlinear  dispersive  material  that  contains  both  frequency-dependent  electric 
conductivity  and  polarization. 

L  INTRODUCTION 

There  has  been  considerable  interest  in  understanding  the  transient  behavior  of  an  ultrafast  laser  pulse  that 
interacts  with  a  nonlinear  dispersive  material.  In  the  last  several  years  many  experimentalists  have  made  use  of 
newly  developed  Kerr  lens  mode-locked  titanium-sapphire  lasers  to  perform  well-controlled  experiments  to  obtain 
accurate  measurements  of  the  transient  behavior  of  ,ultrafast  laser  pulses  in  simple  molecular  liquids  and  solids 
which  are  known  to  exhibit  nonlinear  optical  behavior  [1].  To  better  understand  the  details  of  nonlinear  processes 
that  are  observed  in  the  experiments,  numerical  simulations  have  been  used  extensively  to  reproduce  observed 
nonlinear  effects.  Until  recently  most  computer  simulation  has  been  performed  by  solving  an  approximation  to 
Maxwell’s  equations,  known  as  the  generalized  nonlinear  Schrodinger  (GNLS)  equation  [2],  to  get  information 
about  the  time  evolution  of  the  envelope  of  the  propagating  oscillating  wave  packet  so  that  one  can  obtain  the 
overall  shape  of  the  propagating  optical  pulse.  Because  a  GNLS  equation-based  computer  simulation  does  not 
provide  any  information  about  the  details  of  the  oscillating  waves  inside  the  envelope  of  the  optical  pulse,  there  is  a 
renewed  interest  in  solving  Maxwell’s  equations  directly  without  having  to  rely  on  any  approximations. 

With  the  advent  of  present  day  computers  which  provide  very  fast  execution  times  and  great  quantities  of 
computer  memory,  we  are  at  the  point  where  we  have  enough  computational  power  to  solve  Maxwell’s  equations 
directly  for  nonlinear  dispersive  materials.  Among  recently  investigated  numerical  techniques  that  show  great 
promise  in  achieving  this  goal  is  the  well-known  finite-difference  time-domain  (FDTD)  method  [3].  It  is  based  on 
using  a  simple  differencing  scheme  in  both  time  and  space  to  calculate  the  transient  behavior  of  electromagnetic 
field  quantities.  Because  of  the  simplicity  of  the  FDTD  method,  recent  researchers  have  focused  their  attention  on 
the  numerical  evaluation  of  the  linear  and  nonlinear  convolution  integral  terms  which  appear  in  one  of  Maxwell’s 
equations  (Ampere's  Law).  By  properly  evaluating  these  terms,  many  people  have  successfully  modeled  the 
response  of  linear  and  nonlinear  dispersive  effects  [4-10]. 

In  this  paper  we  consider  the  isotropic  materials  that  exhibit  both  linear  and  nonlinear  polarization  properties, 
specifically  through  the  first-order  (linear)  and  third-order  (nonlinear)  electric  susceptibility  functions,  Xj0(t-r)  and 
Xp3)(t,T,ti,t2),  respectively.  For  such  materials  the  relationship  between  D(t;x)  and  E(t;x)  can  be  expressed  as  [11] 

oo 

D(t;x)  =  e„€„  E(t;x)+e0  £  \E(T;x)X<p,)(t-r)dT 

P  -oo 

x)[ E(t];x)*E(r;x)]Xip)( t,r,tl,t2  )dt2  dtj  dr ,  (1.1) 

p 


1 


where  €„  is  the  electric  permittivity  of  free  space,  is  the  medium  permittivity  at  infinite  frequency,  and  Xp'\t-z) 
is  the  /7th  term  of  the  collection  consisting  of  pmax  time  dependent,  first-order  electric  susceptibility  functions,  where 
/W  is  the  maximum  number  of  terms  which  we  choose  to  consider  for  a  particular  formulation  of  Eq.  (1.1), 
Xp  (t,z,tht2)  is  the  pth  term  of  the  four-time  dependent  third-order  susceptibility  function  which  contributes  to  the 
nonlinear  behavior  of  the  material  and  •  is  the  notation  used  for  the  dot  product  of  vectors.  When  Xp3>(t,z,tht2)  is 
reduced  to  the  single-time  dependent  susceptibility  function,  zj3>(trt2),  by  making  use  of  the  following  Born- 
Oppenheimer  approximation  [11]: 

X<p>(t,T,tI,t2  )=S(T-t,  )S(t-t2  )[%fp>(t2 -tj  )  +  S(t2~t)  d03J] ,  (1.2) 

where  dop  is  a  constant  and  dft)  is  the  Dirac  delta  function,  we  can  show  that  Eq.  (1.1)  reduces  to  an  expression  that 
consists  of  sums  of  convolution  integrals  of  linear  and  nonlinear  terms;  namely, 

oo 

D(t;x)  =  €„€„  E(t;x)+e0  £  j E( T;x)X,pI >( t-t ) dr 

P  -oo 

oo 

+€"  E(t;x)[ E(t;x)*E(t;x) ]^a,3>  +e„  E(t;x)^  J [ E(z;x)»E(z;x)] Z(p3>(t-z)dz  .  (1.3) 

P  P  -oo 

We  also  consider  in  this  paper  the  case  where  the  current  vector,  J(t;x),  which  appears  in  Ampere’s  Law,  is 
represented  by  two  contributions:  the  first  term  is  directly  proportional  to  the  electric  field  vector,  E(t,x),  with  the 
constant  conductivity  coefficient,  6°\  and  the  second  term  is  expressed  in  the  linear  convolution  integral  of  E(t;x) 

and  the  first-order  time  dependent  conductivity  function,  dl;(t).  Hence,  we  express  the  current  vector  in  the  form  as 
shown  below. 


£(t;x)  -  <?<0>  E(t;x)+  ^E(z;x)ain(t-z)dz .  (14) 

Based  on  the  above  expressions,  we  provide  a  general  FDTD  formulation  for  evaluating  the  linear  and  nonlinear 
convolution  integrals  that  appear  in  Ampere's  Law.  We  investigate,  in  particular,  the  case  in  which  the  first-order 
electric  susceptibility  function,  the  third-order  electric  susceptibility  function,  and  the  first-order  time  dependent 
conductivity  function  are  all  expressed  in  the  following  complex  exponential  forms  that  contain  complex  constant 
coefficients: 

X<J>(  t)  —  Re{  aLpexp(—yLpt )}  U(t) ,  (1  5) 

Xip>(t)  =  Re{ aNpLexp(-yNpLt)}U(t),  and  (1  6) 

e<n(t)  =  Re{  pLexP(-dLt)}U(t).  (1  7) 

where  Ref  hs  used^to  represent  the  real  part  of  a  complex  function,  U(t)  is  the  unit  step  function,  and  apL  ,  apNL  , 
P  ’  Yp  ’  Yp  and  6  are  complex  constant  coefficients;  superscripts  L  and  NL  are  used  to  distinguish  between  linear 
an  nonlinear  coefficients.  By  making  the  proper  choices  of  complex  constant  coefficients  and  performing  Fourier 
transforms,  we  can  readily  obtain  the  familiar  Debye  and  Lorentz  forms  of  the  complex  conductivity  and 
permittivity  in  the  frequency  domain. 


II.  FDTD  FORMULATION 


In  light  of  Eqs.  (1.3)  and  (1.4),  we  write  Maxwell’s  equations  inside  the  dispersive  material  as 
d[pH(t;x)  ] 

yt - --Y*E(t;x),  (2.1) 

dD(t;  x)  ...  , 

— yt -  =  Y.  xtL( Ex)  —  <T  0  E{ t;x)—  J_L( t;x),  (2.2) 


with 
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(2.3) 


Dff;jc)=e„e«,  E(t;x)+e0  PLp(t;x ) 

p 

+e„  E(t;x)^jPNpL( t;x)  +  e„  E( t;x)[ E(t;x)»E( t;  x)]^o{03J , 
p  p 

oo 

PLp(t;x)  = 

— OO 

oo 

(f;*)  ■  J7 E(t;x)*E(t:x)]  t-r)dx , 


J_L(t;x)  =  ^E(T;x)CTa)(t-T)dr , 


where  //(X**)  is  the  magnetic  field  vector,  //  is  the  magnetic  permeability,  Pp(t;x)  and  [E( t;x)PpNL( t;x)]  are  related  to 
the  pth  terms  of  linear  and  nonlinear  polarization  field  vectors,  respectively,  and  JL(t;x)  is  the  first-order  time 
dependent  conductivity.  Using  an  FDTD  algorithm,  the  above  equations  can  be  solved  numerically  at  each  time  step 
provided  we  can  handle  J  L(t;x),  Pp(t;x)  and  Ppl(t;x)  numerically.  Therefore,  the  whole  solution  rests  on  the 
question  of  how  to  carry  out  the  numerical  evaluation  of  JL(t;x ),  Pp(t;x)  and  PpNL(t;x)  at  each  successive  time  step. 
For  that  reason,  the  rest  of  this  section  is  devoted  to  the  discussion  of  numerical  formulation  that  treats  JL(t;x)9 
PpL(t;x)  and  PpNL(t;x). 

To  obtain  second-order  accuracy  in  time  in  evaluating  the  convolution  integrals,  E(t;x)  is  approximated  by  a 
piecewise  continuous  function  over  the  entire  temporal  integration  range  so  that  E(t;x)  changes  linearly  with  respect 
to  time  over  a  given  discrete  time  interval  [mAt,  ( m+l)At ],  where  m=0, n,  with  nAt  being  the  current  time  step 
[12,13].  Thus,  the  mathematical  expression  for  E(t;x)  takes  the  following  form  which  can  be  expressed  in  terms  of 
the  electric  field  values,  E}jkm  and  Eykm+\  respectively,  evaluated  at  discrete  time  steps  t-mAt  and  t-(m+l)At  where 
these  two  successive  times  are  obtained  at  the  same, discrete  spatial  location  x=(iAxJAy,kAz)  with  Axf  Ay  and  Az 
being  the  spatial  grid  sizes  in  the  x,  y  and  z  directions,  respectively  (we  use  a  superscript  to  designate  the  discrete 
time  step  and  a  subscript  for  the  discrete  spatial  location): 


—( t*  *  )x~(  iAxJAy.kAz  ) 


=  (K7ik  + 


( En,+I  —Eni  ) 

(t-mAt)  }  for  0  <  mAt  <t  <(  m+1  )At£(  n  +  1  )At ;  and 
At 

for  t<0.  (2 


M.(t;x)xHiAx.jAyjuiz)=  0  for  t<0 .  (2.7) 

For  three-dimensional  FDTD  calculations  in  Cartesian  coordinates,  we  need  to  solve  the  following  discrete  forms 
of  Maxwell's  equations  which  are  obtained  from  Eqs.  (2.1)  and  (2.2)  by  finite  differencing  in  both  time  and  space 
using  the  usual  Yee  algorithm  [3]: 

M(  )x  - m Hff*  )x={^~[(  E?J(km  A.  -(  3m  A  ]  ~[ ( E?(J+m  A  -( E?(H/!)k  A  ]} ,  (2.8) 

M(  A  ~M(  HI*  ),={%[(  E(ni+m  A  ~(E^,m  )z]-~[(  3m  A -( 3m  A  ]}  >  (2-9) 

M(H£*)Z -(i(Hl*  A  ={^[(E?u+m  A  (  E i(j-i/2)k  )x  ]  ^  [  (  E(i+l/2)jk  )y  —tfyjk  )y  ]},  (2.10) 

(  d£'  A  -( D?jk  )x  =  { ~[( H*j%k  A  -( A  *>■  H 

( n+1  )At  (n+l)At 

-<7(0>{  J  ( E(t; x)xMlAx  jAykAz)  A  dt  }-{  j(JL(t;x)x=<i4x  jA>.k/kk)xdt},  (2.11) 
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< D»  -f D •  >*-(%[< i  -< «»-» >•  ] ~i< ";%»» k  -< n<$,»  >,  n 

(n+l)A  (n+l)Ai 

®  (  J  (E(t>X )x=(iAx.jAy.kAz)  )y dt  }  —  {  j"(7  ( I  •  *  h=(  iAx.jAy.kAz )  )ydt  /,  (2-12) 

nAi  nA, 

<  Off'  k-d>i>'‘{±[<  1,  K  -( HZim  HJ&,  >,  »a,  i  11 


( n+I  )At 


( n+I  )At 


-cr,0>{  J  (E(t;x)^iAxjAykAz))zdt  }-{  j(JL(t;x)£=<iAXtjAykAz))zdt}.  (2.13) 

nA>  nAt 

To  evaluate  the  right-hand  side  of  Eqs.  (2.1 1-2.13),  we  first  begin  by  substituting  Eq.  (1.5)  into  Eq.  (2.4)  and 
Eq.  (1.6)  into  Eq.  (2.5)  and  then  differentiating  these  integrals  with  respect  to  time  to  obtain  the  following  first-order 
differential  equations  for  complex  functions  QpL(t;x)  and  QpNL(t;x): 


d&'Jtzx) 

- yt - +YpQp(t;x)  =  aLpE(t;x), 


(2-14) 


dQpL(t;x) 

- -J-t - +Vp  QP  (t;x)  =  aNPL[E(t;x)>E(t;x)].  (2.15) 

From  these  equations,  the  linear  and  nonlinear  polarization  vectors,  Pp(t;x)  and  [E(t;x)PpNL(t;x) ],  can  be  obtained 
simply  by  taking  the  real  parts  of  QpL(t;x)  and  (E(t;x)QpNL(t;x)l  respectively.  Now,  solving  these  two  equations 
exactly  by  using  integrating  factors  exp(ypL  t)  and  exp(ypNL  t),  respectively,  and  then  performing  the  time  integration 
from  nAt  to  (n+I) At,  we  have 

^  nAt+At 

Qp (  n  At  +  At;  x)  =  exp(~yLp  A  t)Qp( nAt; x)  +  aLp  exp( -yLpAt)  J E( t; x ) exp[ -yp( nAt -T )]dt ,  (2.16) 


QpL(  nAt  +  At;  x)  =  exp(  -yNL  At  )QNL(  nAt;  x) 


nAl+At 

+  apLexp(-ypL At)  j[E(T;x)»E(r;x)]exp[-ypL(nAt-r)]dT.  (2.17) 

nAt 

Based  on  the  piecewise  linear  approximation  which  we  have  considered  for  the  temporal  behavior  of  E(t;x),  we  can 
substitute  Eq.  (2.7)  into  the  right-hand  side  (i.e.,  the  inhomogeneous  part)  of  Eqs.  (2.16)  and  (2.17)  for  E(t;x). 
Evaluating  at  spatial  location  x=(iAx,jAy,kAz),  these  equations  result  in  the  following  recursive  relationships  for 
(QP  h/  and  ( Qp  ),/+  in  terms  of  (QpL)iJkn ,  (QpNL)ijkn ,  EiJkn  and  Eijkn+I: 

QLp(nAt  +  At;iAK,jAy,kAzMQLp &*'  =  exp( -  yLp. At ) ( QLp)?jk  +  E^ VLp0)  +  ( E^ -E^  )(VLpl),  (2.18) 

and 

QNpL(nAt  +  At;iAx,jAy,kAz)  =  (QNpL)$'  =  exp( -  yf  At )( QNpL )?jk  +/(£«,  )](¥NpLQ) 

+  2[E"jk.(g^-E^jk)](¥^) 


where 


+[(m'-M!jkwZu;k,-E”Jk)](¥NpL). 


( n+/  W  At  L 

(V , o.o)  =  aP  j  exp[ -■}+(  nAt  +  At -t  )]dr  =  ap  jexp(  —y^T  )  dT—^-~  [  I  —  exp(  ~ypAt )] , 

nAt  0  Y p 

L  At  ^ 

(¥pJ)=~£\(At-T)exP(^LpT)dT  =  ^f{ 1 — r~[l-exp(-yL  At)]}  , 

*  o  Tp  Yp&  P 
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( K Lo > s 0$ 1 J exp( -yf  T )dr=^—-[ l—exp( -y^L At)] , 
o  fp 

NL  At  ml 


(2.22) 


ZlZ  J  r  yv^  v/V^  At  H 

0  I  p  I  p  ** 

NL  At  ml 

(^Z)s~~jj(At-T)2exp(~KLT)dT=^r(1 — T^[1~-yrTlI~exrf~rfLAt) Hi-  (2-24) 

(At)2  l  p  f'p  yNpL At  yfAt 

When  Eqs.  (2.18)  and  (2.19)  are  used  in  Eq.  (2.3)  to  obtain  discrete  forms  of  Eq.  (2.3)  at  two  successive  times 
t-nAt  and  t-(n+])At  and  at  specific  spatial  location  x=(iAx,jAy, kAz),  we  can  express  DiJkn  and  ,  respectively, 
in  the  following  forms: 

ek  =€«€-  %jk+*o  Y<Re((9.%  h*o  gjkJ,Re{(QNpL$jk  }+z„  E]jk[E]jk  .E]Jk  ]%<%>,  (2.25) 


— p7  E];'+e0  J^Refl exp( -yLpAt )[( QLp fa  )  +  E*k ( ¥Lp0 ) + ( Efa1  - E]jk  )(¥LpJ)} 

P 

+  e„  ^Ref  Ey^exp) -ypLAt )( QNpL fa  +  E^I[(Efa  .Efa  )] (¥%  ) 

P 

+2E%'[E“k  .(£%' -E*jk  )]( y/pLi  )+E$'[<e£’ -Efa  ).(£%' -Efa  )]<¥%  )} 

+€„  Efa' IE£‘.e£iYio103J.  (2.26) 

P 

Subtracting  Eq.  (2.25)  from  Eq.  (2.26),  we  obtain  the  following  expression  for  the  left-hand  side  of  Eqs.  (2.1 3-2.13) 
in  terms  of  (QpL)ijkn ,  (QpNL)ijkn ,  Eijkn  and  £,/+/: 

+  G«  ^ Re{ [ exp( -ypAt )-]][( Qp fa  ]  +  E]jk( ¥p0  )+(Efa'  -E*jk  )(¥fa  )} 

P 

+  e0  £ Ref[E];'exp(-y^LAt )  -E*jk](QNpLfa  +g£'[(gfa  *gfa  )](¥%) 

P 

+  2E];‘[E]jk  •(Efa;1 -Efa  )}(¥npl,  )+ Efa1  [(Efa1  -Efa  HE*;' -E*jk )]( ¥pj  )} 
l *■;'[( El .Efa1 )] -Efa  [(Efa  .Efa  )]]%€$.  (2.27) 

P 

To  evaluate  the  second  term  of  the  right-hand  side  of  Eqs.  (2.1 1-2.13),  we  substitute  Eq.  (2.7)  for 
E(t;x)x=(iAxjAy,kAL)  and  obtain  the  following  expression: 

(  n+J  )At 

a<0){  J  E(t;x)x=<iAxJAyJcAz)  dt  }=a<0>  AL(Efa'  +  Efa  )  .  (2.28) 

nAt 

To  evaluate  the  third  term  of  the  right-hand  side  of  Eqs.  (2.1 1-2.13),  we  first  substitute  Eq.  (1.7)  into  Eq.  (2.6) 
and  then  differentiate  the  integral  with  respect  to  time  to  obtain  the  following  first-order  differential  equation  for 
complex  function  JL(t;x): 

+0LlL(t;x)  =  pLE(t;x).  (2.29) 

As  in  the  case  of  linear  and  nonlinear  polarization  vectors,  the  current  vector,  JL(t;x ),  can  be  obtained  simply  by 
taking  the  real  part  of  I_L(t;x)  in  the  above  equation.  Once  again,  we  solve  the  above  differential  equation  exactly  by 
using  integrating  factor  exp(0L  t)  and  then  performing  the  time  integration  from  nAt  to  t  to  get  the  following 
expression  for  complex  function  JL(t;x): 
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(2.30) 


L L(t;x )  =  exp [0L(  nAt -t)]JL(  nAt;x)  +  J. 9Lexp(-0Lt )  j exp(0Lt  )E( T;x)  dr 

nAt 

Again,  we  substitute  Eq.  (2.7)  into  the  above  equation  for  E( t;x)  based  on  the  piecewise  linear  approximation.  Then 
we  integrate  Eq.  (2.30)  from  nAt  to  nAt+At,  evaluated  at  spatial  iocation  x=(iAxJAy,kAz),  to  obtain  the  following 
expression  for  the  third  term  of  the  right-hand  side  of  Eqs.  (2. 11-2.13): 

(n+l)At  { n+l)At 

J—  iAx,jAy,kAz  )dt  =  Ref  J  /  ( t;X  iAxJAv.kAz  )  &  } 

nAt  nA} 


( n+J  )At 

=  Re{ exp( 0LnAt)JL( nAt;iAx,jAy,kAz )  jexp(-0L  t )dt 

nAt 

( n+I  )At  t 

+  PL  J  \exp[0L(r-t)]  E( r;x)x=(iAxJAykAz)  dtdt  } 

nAt  nAt 

(  n+1  )At 

=  Re{  exp(0LnAt )(  JL  )?k  jexp(-0L  t  )dt 


nAt 


( n+1  )At  t 


+  PL  J  J  exp[0L(t-t )][ E"jk  + 


( Ktk  ~Kiik  ) 


nAt  nAt 


At 


( X-nAt )]  dxdt } 


=  Re( J  —  exp( -0LAt )](  JL  )£k  +  E”Jk(tf  )+(£%' -E”jk  (2.31) 


where 


(n+l)Al  t 


\  riTi  T  L 

(£o)  =  PL  J  jexp[0L(T-t)]dTdt  =  -@-^{  1— L-  [l-exp(-0LAt)]}, 


(2.32) 


nAt  nAt 
(  n+I  )Al  t 


a  L  '  n+!  /At  t  l 

=  —  J  ^(T-nAt)exp[0L(T-t)]dTdt=^-AL{L-_L_+ - L - [ l  —  eXp( -0LAt )]} .  (2.33) 

nAt  nAt  0  ^  6  At  ( QL  At  )2 

The  (/  )jjk  expression,  which  appears  in  Eq.  (2.31),  can  be  obtained  in  the  recursive  form  when  we  solve  Eq.  (2.29) 
exactly  for  /  (t;x)  using  integrating  factor  exp(0L  t)  and  then  performing  the  time  integration  from  t-(n-l)At  to 
t=nAt  while  making  use  of  the  piecewise  linear  approximation  for  E(t;x)  to  evaluate  the  integral,  which  arises  from 
the  right-hand  side  of  Eq.  (2.29).  The  result  is  that  we  obtain  the  following  recursive  relationship  for  (JL) ijk  which  is 
expressed  in  terms  of  the  previous  time  step  values  of  (/%*"' ' ,  Eijk''  and  Eijk : 

=exp(-eLAt)aL)ik' +*£(&+(&  -&)<#)}, 

where 

nAt  At  . 

P  r 


(Co)-PL  J  exp  [~0L(  nAt-T  )]  dr-  pL  jexp(  -0Lr  )dr=^-[  ]-exp(  -0L  At )] 

(n-l)At  0  & 

/?  ^  ^  ^ 

j(At-T)exp(-0LT)dT  =  ~f  1 - \—[l-exp(-0lAt)]}  . 


(2.34) 

(2.35) 


(2.36) 


Finally,  when  Eqs.  (2.27),  (2.28)  and  (2.31)  are  substituted  into  Eqs.  (2.1 1-2.13),  we  obtain  the  following  three 
coupled  nonlinear  cubic  equations  which  we  need  to  solve  for  Eijkn+I  [i.e.,  (EiJkn+,)x ,  (EiJkn+,)y  and  (Eijkn+I),  ]  in  terms 

of  known  quantities  (Qp  )yk  ,  (QpNL)ijicn  ,  (X')itk  ,  £</*"  and  tLjk"+/'  which  are  calculated  at  previous  time  steps  t=nAt 
and  t=(n+%)At: 

ao*  +aP  E$‘  )x  +(E!jZ'  )J a2x(Eg'  )x  +a2y( Eyk'  ),  +  a2z(E%'  )z  } 

)J  [(£?;’  )x  7 2  +[(E£<  )y  ] 2  +[(£$'  )z  ]2  }  =  0,  (2.37) 
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(2.38) 


+a,(E?kl ),  +(E%'  ),{ a2x(Eg'  )x  +a2y(E )y  +a2l(E$'  )z  } 

+  a3(E£>  )J  [(E?;‘  )x  ] 2  +[(£$’ ),  ]2  +[(Eyk'  )z  ] 2  }  =  0, 
a0z  +a,(E I?  )z  HE?;1  )J a2x(E?;‘  )x  +a2y(E £'  )y  +a2z(E?;’  )z  } 

+h(e;'  )J  he?;'  )x  i2 +[<£?;'  )y  p  +[(£?;'  )z  p}=o, 

where  a0x ,  a^. ,  ,  a, ,  a2x ,  a2y ,  a2z  and  a3  are  given  by 

aoH(^[(  h?;»,a)  )y  h  h?;;,/2)  a  7  ~U  n$U  a  -( n;%k  a  77 + ^(0J  *5* 

W  )?jk  ]x  }+( E?jk  )x  Refftf  )-(&)) 

0 

-  (££*  )x+€„  ^Re{ [exp(-  ylpAt )- 1 ][( QLp%k  ]x  }+e0  (E?Jk  )x^Re{(¥Lpj0  )-(yLpJ )} 

P  P 

(E?Jk  )^Re{«fpL;  }-e„  (  E?jk  )J[(E”jk  )x]2+[(  E?jk  )y]2+[(E?jk  )zf  <$, 

P  P 

a0\  ~ (  ^  j+i/i) jk  )Z  -(  H(  i-V2)}k  A  7 -%[ (»%£»  A  -( HSSHw  A  ]}+cr(l>>~( Eyk  A 

+  Re{  Jr[J-exp(-eLAt)J[(jL  fa  ]y  }+(Enjk  )y  Re{( )-( )} 

0 

-  (E!‘k  k+B(>  X Ref  [exp(-yLpAt  )-l  l[(Q%ly  }+e„  (E?jk  )y^Ref(¥p.o  )-(¥pj  )} 

P  P 

-e„  ( E?Jk  )^Re{(QNpLf:jk}-B0  ( E?Jk  )y[[(E?jk  )xf+[(E»k  )y]2+[(E?jk  )z]2  c$ , 

P  P 

aozM~[(  A  -  (  h?<hw  u~[(  Hfiwjt  A  -( Hu-i)jk  A  77 + y  (  E?jk  A 

+  Ref  l[I-exp(-eLAt)][(jL  fa  ]z  }+(E?jk  )z  Re{(4L0)-(tf)} 

& 

-  G„e„(E?jk  )z+et>  ^Ref  [exp(-yLpAt)-l][(QLpfa  ]z  }+e„  ( E“k  )z  ^Ref(¥^0H¥pj  )} 

P  P 

~e0  (  E?jk  A  ^Re{(QNpL]?jk  }-e„  (E?jk  A  [l(E?jk  )x]2+[(E?jk  )y]2+[(E?Jk  )zf  «£>  , 


(2.41) 


(2.42) 


a,  Se„e^  +(T«»*L+Re{(ti)}+€o  X /?*/(<;  )}+e„  J^Re, { exp(  -  yNpL  At )( QNpL  fa  } 


P  P 

+  €„  [ [(E?jk  )x  ]2  +[( Eyk  )y  ]2  +  [( Ejjk  ).  ]2  J^RefOr^o  )-2(¥NpL,  )+(¥NpL2  )} ,  (2.43) 

P 

a2x= 2eJE?jk  )xy£Ref(^L,)H¥pL2)L  (2-44) 

P 

H  (E?jk  A-X  ReUvf,  )H¥pL2  )} .  (2-45) 

P 

<*2z  s2e„  (  E?jk  AX  Ref  ( ¥pj )~( ¥pL2  )} ,  (2-42) 

P 

a3  se„  ^Re{(¥pL2  )}+ X  •  (2-43> 

P  P 


Eqs.  (2.37),  (2.38)  and  (2.39)  can  be  solved,  respectively,  for  (Eijkn+I)x  ,  (EiJkn+I)y  and  (EiJkn+1)z  by  using  any  standard 
root-finding  numerical  technique  for  a  set  of  nonlinear  equations  [14].  One  technique  that  is  suitable  for  our  problem 
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is  the  iterative  nonlinear  Newton-Raphson  method  using  (Eijkn)x  ,  (Eijkn  ),  and  (EiJkn  )z  as  the  initial  guess  for  the  start 
or  the  iterative  procedure. 

In  the  above  formulation,  we  only  need  to  consider  updating  Eqs.  (2.8-2. 1 1),  (2.18),  (2.19),  (2.34)  and 
(2.37-2.39),  respectively,  for  ,  (QpL)ijkn+l ,  (QpNL)jjkn+l ,  (JL)ijk"  and  Eijk"+ 1  at  each  time  step  in  order  to  carry  out 
a  complete  three-dimensional  computer  simulation  of  the  electric  field  response  in  nonlinear  dispersive  materials 
that  contain  both  frequency-dependent  electric  conductivity  and  polarization  terms. 

For  the  purely  linear  dispersive  case,  a2x ,  a2> ,  a2z  and  a3 ,  as  well  as  some  terms  appearing  in  a0x ,  a0y ,  a0z  and  a, , 
turn  out  to  be  zero.  In  this  case  we  can  solve  for  Eijkn+I  directly  without  having  to  rely  on  the  numerical  root  finding 
technique  as  discussed  in  greater  detail  in  the  published  literature  [15-20]. 


III.  CONCLUSIONS 


Based  on  the  FDTD  approach  we  presented  here,  we  can  solve  Maxwell’s  equations  directly  for  the  propagation 
ot  electromagnetic  waves  in  linear  and  nonlinear  dispersive  media  that  exhibit  the  frequency-dependent  electric 
conductivity  and  polarization.  Because  of  the  piecewise  linear  approximation  we  used  for  the  time  dependent  part  of 
the  electric  field  vector,  our  approach  provides  second  order  accuracy  in  time.  In  addition,  our  approach  retains  all 
the  advantages;  of  the  usual  first-order  discrete  recursive  convolution  approach,  such  as  fast  computational  speed  and 
efticient  use  of  the  computer  memory. 

We  .J’a.ve  shown  that  il  is  critical  t0  express  the  first-order  (linear)  conductivity,  the  first-order  (linear) 
susceptibility,  and  the  third-order  (nonlinear)  susceptibility  functions  in  the  exponential  forms  in  the  time  domain  in 
order  toobtain  the  recursive  feature  in  our  FDTD  algorithm.  The  most  important  result  we  have  shown  in  this  paper 
is  that  FDTD  formulation  for  nonlinear  dispersive  materials  results  in  having  to  solve  three  coupled  nonlinear  cubic 
equations  for  the  three  components  of  the  electric  field  vector  at  each  time  step  as  compared  to  just  solving  the  linear 
equations  in  the  case  of  linear  dispersive  materials. 
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